Simulation of Binary Black Hole Spacetimes with a Harmonic Evolution Scheme 
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A numerical solution scheme for the Einstein field equations based on generalized harmonic co- 
ordinates is described, focusing on details not provided before in the literature and that are of 
particular relevance to the binary black hole problem. This includes demonstrations of the effec- 
tiveness of constraint damping, and how the time slicing can be controlled through the use of a 
source function evolution equation. In addition, some results from an ongoing study of binary black 
hole coalescence, where the black holes are formed via scalar field collapse, are shown. Scalar fields 
offer a convenient route to exploring certain aspects of black hole interactions, and one interesting, 
though tentative suggestion from this early study is that behavior reminiscent of "zoom-whirl" or- 
bits in particle trajectories is also present in the merger of equal mass, non-spinning binaries, with 
appropriately fine-tuned initial conditions. 



I. INTRODUCTION 

The past year has seen a remarkable leap in progress 
toward a numerical solution of the binary black hole 
problem. The first stable evolution of an entire merger 
process — orbit, merger, ringdown and gravitational wave 
extraction — was presented in Q using a harmonic for- 
mulation of the field equations. Shortly afterwards two 
groups independently achieved similar success using 
a modified form of the BSSN (or NOK) 0,1111] formula- 
tion of the field equations. These first results all focused 
on the merger of equal mass, non-rotating black holes, 
and follow-up studies 7, 8] are now honing in on a con- 
sistent picture of the gravitational waves emitted by such 
an event. Recently, similar techniques to those employed 
in 0, 13 were successfully used to study unequal mass 
mergers H and provide estimates of the "kick" ve- 
locity imparted to the final black hole, and in the ef- 
fects of black hole spin, aligned and anti-aligned with the 
orbital angular momentum, were studied, demonstrating 
that in the aligned case some orbital "hang-up" occurs 
to radiate away the excess angular momentum. 

One reason why the new results seem like such a leap 
is that progress during earlier years had been rather slow 
and arduous. For example, it had taken roughly 7 years 
from the first simulation of a fraction of a non head-on 
collision to a full orbit J 3], with many groups making 
advances along the way. Much of the effort was also 
(and still is) focused on understanding the underlying 
structure of the field equations of general relativity, and 
why they are so problematic to discretize successfully in 
many cases. Nevertheless, the manner in which advances 
had been made over the pass decade suggested a similar 
series of "baby steps" toward a future solution, which is 
why the results of the past year have been so exciting. 

The primary purpose of this paper is to describe in 
detail certain aspects of the generalized harmonic (GH) 
evolution scheme that arc of relevance to the simulation 
of binary black hole (BBH) spacetimes, as first presented 
in P, Il4| . This is not a comprehensive overview of the 
code or technique, and could be viewed as a supplement 



to 0, Q (and see also 0| for an excellent description 
of GH evolution). A secondary goal is to relate progress 
on an ongoing effort to understand the nature of BBH 
spacetimes produced by scalar- field collapse^. It is ar- 
guable how relevant such spacetimes may eventually be 
to gravitational wave detection efforts, nevertheless they 
offer an easy route to explore a wide range of BBH pa- 
rameter space. One interesting though tentative result 
from this early study is that the zoom-whirl type be- 
havior seen in geodesies around black holes may also be 
present in BBH orbits involving comparable mass com- 
ponents. In particular, what is shown is that for equal 
mass binaries on non-circular orbits, and with appropri- 
ate fine-tuning of the initial conditions, the black holes 
approach one another, "whirl" around for several orbits, 
then separate again. 

The outline of the rest of the paper is as follows. In 
Sec.^a brief overview of the Einstein field equations in 
harmonic form with dynamically evolved source Junctions 
and constraint damping terms is given. Sec. IIIII contains 
a description of parts of the numerical code that seem 
to be important for stable evolution of BBH spacetimes 
with this scheme, including excision and numerical dissi- 
pation. The effects of constraint damping and the choice 
of source function evolution equations on a select set of 
simulations is presented in Sec. IIVI Results from the 
inspiral and merger of scalar field generated BBH space- 
times are given in Sec.0 followed by concluding remarks 
in Sec. EH 



II. GENERALIZED HARMONIC 
COORDINATES 

In this section a brief description of generalized har- 
monic coordinates, in the specific form that is discretized 



^ a study of Cook-Pfeiffer quasi-circular inspiral initial data ll6l 
evolved using the g ener alized harmonic evolution scheme will be 
presented elsewhere 1 17l 



2 



in the code, is given. For a thorough descri pti on and the 
history of the use of these coordinates see [13 • 

Consider a spacetime describe by the hne element ds, 
metric tensor gab and coordinates x° 

ds^ — gabdx°'dx^ , (1) 

and the Einstein field equations in the form 

Rab = StT (^ab - \9abT^ , (2) 

where Rab is the Ricci tensor, gat is the metric tensor, 
Tab is the stress energy tensor with trace T, and units 
have been chosen so that Newton's constant G and the 
speed of light c are equal to 1. The Ricci tensor is defined 
in terms of the Christoffel symbols FJ^^ 

Kb = ^ff'"' [9ae,b + gbe,a " 9ab,e] (3) 

via 

p j^d T^d I T^e j^d T^e j^d f a\ 

^ab — i ab,d ~ db,a ^ ^ ab^ ed ~ ^ db^ ea \^} 

The notation f^a and daf is used interchangeably to de- 
note ordinary differentiation of some quantity / with re- 
spect to the coordinate x". 

Harmonic coordinates are defined by the following set 
of four conditions on the four spacetime coordinates: 

□x" = 0, (5) 

where □ is the usual covariant scalar wave operator: 

= VbV^x- - ■ (6) 

Generalized harmonic ( GH) coordinates introduce a set 
of arbitrary source functions into the definition (jSJ : 

Ux' = H". (7) 

Note that in contrast to harmonic coordinates (0, GH 
coordinates Q are not conditions on the spacetime 
coordinates — any geometry in any coordinate system can 
be expressed in GH form with defining the corre- 
sponding source functions. 

Using the definitions above, and defining He — gacH"', 
the Einstein equations ||2Jl can be rewritten in the follow- 
ing equivalent form, which shall be referred to as the gen- 
eralized harmonic decomposition of the field equations: 

■^9'"^9ab,cd + g'"^(,a9b)d,c + H(a,b) - Hd^tb 

+rgdFi = -Stt [t^p ~ \9c.pT^ (8) 

The utility of the GH decomposition in a numerical 
evolution (as first carried out in [TsL HqI l20l |. and 



also 0, 121I ) comes from considering the He as in- 
dependent functions] ||SJ| then becomes a set of ten man- 
ifestly hyperbolic equations for the ten metric elements 
gab- As He are now four independent functions, one needs 
to provide four additional, independent differential equa- 
tions to solve for them, schematically written as 

CcHc — (no summation). (9) 

Cc is a differential operator that in general depends on 
the spacetime coordinates, metric and source functions. 
To complete the specification of the system one needs 
to provide evolution equations for matter sources, which 
couple to the field equations through the stress energy 
tensor. The only matter field considered here is a mass- 
less scalar field $ that satisfies the wave equation 

□$ = 0, (10) 

and has a stress energy tensor 

Tah - 2$,a$^b - .gafc*x$'', (H) 

A solution to JHJ, © and (|l()|l will be a solution to the 
Einstein equations as long as the GH condition (0) is 
satisfied for all time. In theory, this is straight-forward 
to achieve. Define the GH constraint functions C° by 

C" = H" - Ux". (12) 

Any solution to the Einstein equations must have C" = 0. 
Using the contracted Bianchi identity and conservation of 
stress energy, one can show that C" satisfies the following 
homogeneous wave equation 

UC" = -R^bC". (13) 

Therefore, what needs to be done to ensure that ©, © 
and (|10|l satisfy the Einstein equations for all time t is to 
choose initial conditions for gab and Ha such that C"" — 
and dtC^ — aX t — 0, boundary conditions on gab and 
Ha such that C° = on the boundary for all time, and 
couple this to matter that conserves stress energy. Then 
(|13|l guarantees that C° = throughout the interior 4- 
volume of the spacetime. 

Of course, things are not as simple as this in practice. 
In a numerical simulation one can only satisfy the con- 
ditions described in the preceding paragraph to within 
the truncation error of the numerical scheme. Again, 
in principle this is not a problem, however it turns out 
that in many situations the truncation errors grow too 
rapidly to achieve useful results given limited resolution. 
There is good evidence that one of the reasons for the 
rapid growth of truncation error is not a poor choice of a 
numerical algorithm, rather H13I) admits rapidly growing 
solutions (so called "constraint violating modes") given 
initial data where C° is of order the truncation error. An 
effective way of dealing with this problem is the addition 
of constraint damping terms to the equations. 
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A. Constraint Damping 



for binary black hole evolutions, namely: 



As introduced by Gundlach et al.^], building on the 
idea of so-called A systems |2J], the field equations with 
constraint damping are 



cd , cd I LJ TJ T'd 

zg gab,cd + g {,agb)d,c + J^{a,b) - -tidl- at 



<1 



1 



= -St: ^T^p - -g^pT j , (14) 

where k is a parameter multiplying the new terms, Ua a 
timelike vector, and Ca are the constraints Since 
the new terms are proportional to the constraints, a solu- 
tion to the Einstein equations (0) will also be a solution 
to H14|) : furthermore, for general solutions of H14|) the 
constraints still satisfy a homogeneous wave equation 



^{bfja) 



(15) 



and thus the prescription outlined in the previous section 
for obtaining valid solutions to the Einstein equations can 
also be applied using 114() instead of ((SJ. Gundlach et 
al. have shown that for perturbations about Minkowski 
spacetime, all finite wavelength solutions to (|15|l are ex- 
ponentially damped in time. It is not known whether 
similar modifications to other forms of the field equations 
can be made that will also have this desirable constraint 
damping property, though in it was shown how to 
apply constraint damping to the ZA formalism |25l |. and 
recently Lindblom et al.jlSi] described a first order sym- 
metric hyperbolic version of the GH decomposition with 
constraint damping. In a first order form of the equations 
additional constraints are introduced that could also ex- 
hibit poor behavior with regards to being satisfied in 
a numerical evolution, though |T5l | describe an effective 
method for dealing with this problem. 

In |22j it was suggested that in (|14|) could be any 
timelike vector field; here, for simplicity, n° is chosen to 
be the unit timelike vector normal io t = = const, 
surfaces. Specifically, Ua = —adat, with a — \/-~^Jg^ 
being the usual lapse function. 



B. Source Function Evolution 

Within the GH decomposition one can think of the 
source functions Ha as representing the four coordinate 
degrees of freedom available in general relativity. There 
are many conceivable ways of choosing Ha — see \\M f or a 
more general discussion of some possibilities, and |26ll2^ 
[28, 29, 30, 31, 32] for coordinate conditions that might 
be readily applicable to GH evolution. The focus here 
will be on one class of equations that have proven useful 



a'' 



^2Ht.yn\ H, = 0. (16) 



This equation for Ht is a damped wave equation with 
forcing function designed to prevent the lapse from de- 
viating too far from its Minkowski value of 1. The pa- 
rameter ^2 controls the damping term, and ^2;'7 regu- 
late the forcing term. The GH numerical code descrided 
here tends to develop instabilities in pure harmonic gauge 
when the lapse function becomes on the order of 0.1 near 
the horizon of black holes; equation H16(l prevents this 
from happening. It is unclear whether the instabilities 
are numerical in nature or an indication that the har- 
monic gauge is developing a coordinate pathology. How- 
ever, this is not a crucial question to answer at the mo- 
ment given the limited availability of computational re- 
sources to fully explore the issue, and that Ijltill works 
well. In Sec. IIVI examples of the effect of (|16|) with a few 
different parameters are shown. 



C. Boosted Scalar Field Initial Data 

An easy way to produce binary black hole "initial" 
data is to use scalar field gravitational collapse. At i = 
one begins with two Lorentz boosted scalar field profiles 
with initial amplitude, separation and boost parameters 
to approximate the kind of orbit that the black holes, 
which form as the scalar field collapses, will have. On 
roughly the light crossing time scale of the orbit the rem- 
nant scalar field that did not fall into either black hole 
propagates away from the vicinity of the orbit, leaving 
behind a good approximation to a vacuum BBH system. 

The procedure used to calculate the initial geometry 
is based on standard techniques '33*1 . One starts with a 
metric written in ADM 34] form 

ds^ = -a^dt^ + h,j {dx' + (3'dt) {dx^ + (3^dt) , (17) 



where as before a is the lapse, (3^ the shift vector and 

const, surfaces. The 
is defined as 



hij is the metric intrinsic to t 
extrinsic curvature of hi 



-h/h.r^^ni. 



(18) 



In the ADM decomposition the Hamiltonian and momen- 
tum constraint equations take on the following form: 



VbKa'' ~VaK ^SirJa, 



(19) 
(20) 



where '^^i? is the trace of the Ricci tensor of the spatial 
metric, K is the trace of the extrinsic curvature, p is 
the energy density and Ja the momentum density of the 
matter in the spacetime: 



p = Tabu^n'' 



(21) 
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Ja = -Timh\n'^. (22) 

For the results described here the following initial condi- 
tions for the geometry are chosen dX t — 0: the spatial 
metric and its first time derivative is conformally flat: 

hij\t=o ^ ■ Vi] (23) 
dthij\t=o ^ dt-ip ■ rjij (24) 

where V' is the conformal factor, and rjij is the flat Eu- 
clidean metric; the initial slice is maximal: 

K = 0, dtK^Q (25) 

and the initial slice is harmonic: 

Ht = 0, H, = 0. (26) 



Minkowski rest frame ds^ = —dt'"^ + dx'"^ + dy'"^ + dz'^: 

<fit' ^0,x',y',z') = Aexp(^-^^ 

dt'Ht' ^0,x',y',z') =. 0, (28) 

where A and A are constant parameters, and r' — 
\J x'"^ + y'^ + Then, a Lorentz boost with velocity 
V is applied in the direction u to this pulse, mapping 
{t',x',y',z') = (0,0,0,0) to {t,x,y,z) = {O,xo,yo, zq), 
where (xq, yg, Zq) is the initial location of the pulse in sim- 
ulation coordinates. Note that this initial data scheme 
can only produce a black hole (when the amplitude A is 
sufhciently large) with approximately the velocity v, as 
no back reaction effects are taken into account^. 



When the conditions (|23I25|) are substituted into the 
constraints H19I20|I . the Hamiltonian constraint becomes 
an elliptic equation for the conformal factor ?/;, and the 
momentum constraints become elliptic equations for the 
components of the shift vector When the maximal 
slicing condition (|25|l is expanded in terms of the metric 
via its definition (|18|l . and using the above conditions, 
an elliptic equation for the lapse a results. This set of 
five coupled elliptic equations is solved using multigrid 
techniques as described in [l^ . For boundary conditions 
the metric is assumed to be the Minkowski metric (and 
these conditions are applied exactly, as a coordinate sys- 
tem compactified to spatial infinity is used). For scalar 
field collapse initial data no "inner" boundary conditions 
are needed as there are no black holes present in the ini- 
tial slice. After the elliptic equations are solved for a, ip 
and H25() is algebraically solved for the initial value 
of dttjj, and similarly the conditions l|26() together with 
the definition ((JJ) is used to solve for the initial values of 
dta and 9(/3'. Once all the initial values and first time 
derivatives of the metric in ADM form (|17|l are known, it 
is straight-forward to convert these to initial conditions 
for the 4-dimensional metric gab needed for the GH sys- 
tem of equations p4|l . In addition to (|26|l . the first time 
derivative of Ht is needed for initial conditions to H16|l . 
and this is chosen to be: 

dtHt - 0. (27) 

Note that the above procedure produces initial data that 
is consistent with the GH constraints H12|l — 0, and 
dtC"' = at t = 0: C° = explicitly as this definition is 
used to provide the initial time derivatives of the lapse 
and shift, and initial data that satisfies the ADM con- 
straints (|19I20() will, to within truncation error, satisfy 
dtC =0 then 35] . 



III. NUMERICAL SIMULATION WITH 
GENERALIZED HARMONIC COORDINATES 

Here a very brief overview of the numerical code imple- 
menting the GH system of equations (|10|1 . (|14|l and 
is given, focusing on a few technical details and miscella- 
neous information related to binary black hole simulation 
that are not discussed in The numerical code has 

the following features: 

• Equations l(Tn)l. ifHI) and ((TE|l are discretized using 
standard second order accurate finite difference stencils. 
The evolved variables are the ten covariant metric ele- 
ments gab, the four source functions Ha, and the scalar 
field A three time level scheme is used, where un- 
known quantities at the most advanced time level are 
solved for using a pointwise Newton-Gauss-Seidel relax- 
ation, given the known quantities at the two past time 
levels. 

• The constraint equations H19|l and H20|) are solved 
with a Full Approximation Storage (FAS) adaptive multi- 
grid algorithm. 

• Excision is used to eliminate the singularities inside 
black holes. The excision surface is an ellipse whose shape 
is a shrunken version of a best-fit match to the shape of 
the apparent horizon (AH) of the black hole. The ellipse 
is shrunk, by typically 0.5 to 0.9, so that the excision 
surface is always some distance inside the AH. The ex- 
cision surface is redefined each time the AH is search 
for, which is frequently enough that the excision surface 
never moves by more than a single grid point between 
searches. This is to ensure stability as extrapolation is 
used to initialize newly "repopulated" grid points, and 
extrapolation tends to be unstable if more than a single 
point inward from the excision surface is repopulated at 
any one time. 



The scalar field initial data is constructed by first tak- 
ing a spherical, time symmetric Gaussian profile in a 



^ for a boost parameter of order 0.2, as used in the simulations 
described in Sec.|y] the estimated velocities of the resultant black 
holes can differ from v by as much as 20 — 30%. 
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• A spatially compactified coordinate system is used 
so that exact Minkowski spacetime boundary conditions 
can be placed on all variables. 

• A combination of adaptive and fixed mesh refine- 
ment is used to efficiently resolve the relevant length 
scales in a given simulation. The grid hierarchy is adap- 
tive in the "near-zone" to track the motion of the black 
holes through the domain, while outside of this in the 
"wave-zone" the mesh structure is kept fixed. This ap- 
proach, rather than full adaptivity, is for computational 
efficiency — the simulations would take a prohibitively 
long time if the mesh structure were allowed to track 
the outgoing wave train. The software libraries used to 
implement the parallel adaptive infrastructure are pub- 
licly available |36l |. though at present the documentation 
is sparse. 

• Numerical dissipation is used to control high fre- 
quency "noise" that often arises in such adaptive simula- 
tions. Furthermore, dissipation is necessary to eliminate 
a high frequency instability that would otherwise occur 
near black holes, where the local lightcones start to tip 
in towards them |33 . 

The remainder of this section contains a discussion of 
how certain properties of the solution are measured, and 
a few technical details of the code: position dependent 
dissipation and constraint damping parameters, and how 
some robustness problems in the apparent horizon finder 
are dealt with. 

A. Measurement of Solution Properties 



nary the Newman-Penrose scalar ^'4 is used, with the null 
tetrad constructed from the unit timelike normal n'^, a 
radial unit spacelike vector normal to r = constant coor- 
dinate spheres, and two additional unit spacelike vectors 
orthogonal to the radial vector'^. Far from the source, the 
real and imaginary components of ^'4 are proportional to 
the second time derivatives of the two polarizations of the 
emitted gravitational waves. 

The following formula is used to estimate the total 
energy E emitted in gravitational waves, 

-dT^A^ J ^ " y ^^'^^ ■ J ^^^^ 

where 'J'4 is the complex conjugate of ^'4, and the surface 
integrated over in 1)31(1 is a sphere of constant coordinate 
radius R. However, before applying this formula ^'4 is 
filtered by eliminating all but the £ = 2, m = ±2 spin 
weight —2 spherical harmonic components -2^,m(^j0)i 
which are the dominant modes of the gravitational 
wave^. This eliminates some high-frequency "noise" that 
is present in the bare waveform (primarily from mesh- 
refinement effects), however the integrated energies in 
the filtered vs. unfiltered waveform do not differ by more 
than 5% at most in a typical simulation. The plots of 
waveforms given in Sec. IVIshow the dominant harmonic 
components of '^4 as a function of time, calculated over 
a given coordinate sphere of radius r: 

_2Q,,„(i, r)= f ^i{r, t, 9, 0) • _2r,,,„(0, 0) dn (32) 



Here a brief summary is given of how black hole prop- 
erties are measured and gravitational waves are extracted 
from a numerical solution. 

At present only the apparent horizons and not event 
horizons of black holes are searched for in the solution. 
Black hole masses are estimated from the AH area A and 
angular momentum J, and applying the Smarr formula: 

M = ^Ml + ,P/{AMl), M„ = ^JA/I&n. (29) 

The angular momentum of horizons are calculated us- 
ing two methods. First, by using the dynamical horizon 
framework [s^ , though assuming that the rotation axis 
of the black hole is orthogonal to the z = orbital plane, 
and that each closed orbit of the azimuthal vector field 
lies in a z = constant surface of the simulation. Due 
to the symmetry of the initial data, these assumptions 
are probably valid, though this will eventuall y n eed to 
be confirmed. The second method, following j^^f) is to 
measure the ratio Cr of the polar to equatorial proper 
radius of the horizon, and use the formula that closely 
approximates the function for Kerr black holes: 

a w \Jl- (2.55a - 1-55)^ (30) 
To calculate the gravitational waves emitted by the bi- 



B. Position Dependent Dissipation and Constraint 
Damping 

With this evolution scheme numerical dissipation is 
essential to control what would otherwise be a high- 
frequency instability near black holes — see [2lll37l |. Also, 
dissipation is useful in suppressing spurious high fre- 
quency components of the numerical solution that are 
sometimes generated at mesh refinement boundaries. 
The Kreiss-Oliger style dissipation employed (see T^) 
converges away in the continuum limit, though with typ- 
ical resolutions used in these simulations the dissipation 
can cause noticeable degradation in gravitational waves 
measured far from the source. Experiments suggest a 
dissipation parameter e (0 < e < 1) of at least 0.3 — 0.5 is 
needed for stability near black holes, though far from the 
black holes much less is needed to control mesh refine- 
ment "noise" . Therefore, in these simulations a position 



^ At this stage all the subtleties in cho osing an "appropriate" 

tetrad are ignored — see for example 1401 
^ note that in the waveform shown was unfiltered, and for the 

energy calculation the filtering was done with scalar spherical 

harmonics 
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dependent dissipation parameter is used which is as large 
as needed near black holes for stability, and then drops to 
a smaller value in the wave zone for improved accuracy of 
the gravitational waveform^. Specifically, inside an AH 
a value of e = 0.5 is used on the excision surface (which 
is between 50% and 90% the size of the AH), then e is 
decreased to 0.35 linearly with coordinate distance from 
the center of the AH to its surface. Outside any AH e 
is set to 0.35 at r = (r = \J ^ ^ is coordinate 
distance from the origin), then linearly interpolated to 
0.05 at r R:! lOMg (A/q is the sum of initial black hole 
masses); for r greater than this e is kept fixed at 0.05. 

In certain situations it may be necessary to use a po- 
sition dependent constraint damping parameter, i.e. re- 
define K in (|14() to be K{x,y,z). As demonstrated in 
Sec. lIVI K needs to be of order 1/L, where L is the small- 
est lengthscale present, to be effective. However, if k is 
much greater than 1 numerical instabilities may develop 
near the outer boundaries of the domain. In those sit- 
uations the following function for k{x, y, z) can cure the 
problem: 

y, z) = [(1 -t- x^){l + y2)(i + ^2)] -m/2 ^ ^33^ 

where m is a positive integer and is a positive constant. 
For all results presented here m = 0. 



C. Excision and Apparent Horizons 

Black hole excision is used to deal with the physical sin- 
gularities that occur inside of black holes. Excision is the 
placement of an artificial boundary around the singulari- 
ties, though inside each black hole. It is possible to do so 
because causality will prevent these unphysical interior 
boundaries from affecting the exterior solution. In the 
code described here all characteristics of the differential 
equations being solved are assumed to be directed into 
the boundary. Thus, no boundary conditions are placed 
on the excision surface; rather, the difference equations 
are solved there using modified finite difference operators 
that do not sample points in the computational domain 
that are excised (see for details). 

The surface that defines the excision boundary is 
guided by the apparent horizon as described at the be- 
ginning of Sec. mil An AH is search for with sufficient 
frequency that the black hole does not move by more than 
a single grid point in any direction between AH searches. 
A flow method is used to find the AH, using the shape of 
the previously found AH as an initial condition. In the 
flow method a tolerance tq is chosen, and the flow equa- 
tion is iterated until jflj^j ^ tq, where \0\t^ is the ^2 norm 
of the expansion of the AH. This works well during most 



^ the simulations shown in used a constant value of e throughout 
the domain 



of the evolution, though near the time of the merger the 
flow method tends to become unstable. If an AH is "lost" 
the simulation crashes shortly afterwards as the excision 
surface is unable to track the motion of the black hole, 
and soon a domain exterior to the black hole is excised. 
One of the factors that seems to cause this behavior in 
the AH finder is relatively poor resolution of the under- 
lying numerical solution. So one cure would be to allow 
additional refinement about the black holes. However, 
experiments have shown that increasing the resolution 
only near the black holes does not increase the accuracy 
in the overall solution, so this would be a computation- 
ally expensive solution merely to help the AH finder. Of 
course, a better solution would be to develop a more ro- 
bust AH finder (see 42, "43] for a review of methods), and 
this path will be pursued in future work. 

In leu of a more robust AH finder the following "tricks" 
are used to push the evolution through the merger point. 
When the fiow method becomes unstable, what typi- 
cally happens is \9\t^ decreases to some minimum value 
Tm > To, then begins to increase and eventually diverge. 
If T„i is not too much larger than tq^, the surface corre- 
sponding to = Tm is used and evolution continues. 
On occasion just before a merger Tm does become too 
large; in that case the motion of the AHs are extrapo- 
lated until an encompassing AH is found, using previ- 
ously measured angular and radial velocities of the AHs. 
This extrapolation is rarely need for more than w 5Mo. 

Note that when the AH finder fails and either a less 
accurate or extrapolated shape is used, this shape only 
guides the position and orientation of the excision sur- 
face, not its size. This is important for stability and 
to prevent the AH robustness problems from adversely 
affecting the exterior solution, as typically the approxi- 
mate AH grows with time, possibly even moving outside 
the event horizon. However, as the AH shape is used to 
measure the mass and angular momentum of the black 
holes the corresponding plots (see Sec. II VI and Sec. 01 do 
refiect the problems of the AH finder. 



IV. CASE STUDIES OF THE EFFECT OF 
CONSTRAINT DAMPING AND CHOICE OF 
GAUGE 

In this section some results are given on the effect of 
constraint damping and choice of source function evolu- 
tion in a dynamical simulation. Little is known about 
these two topics in general, namely, whether constraint 
damping will work in generic 3 + 1 evolutions to sup- 
press the growth of constraint violating modes, or what 
classes of source function evolution equations could be 



typically if Tm, < 50ro , where the factor of 50 was chosen from ex- 
periments in "normal" situations showing that the corresponding 
AH shape differs from the actual AH shape by at most around 
20% in size 
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used to achieve various well-behaved slicings of dynami- 
cally evolved spacetimes. These questions certainly can- 
not be answered with a few cases studies, and that is not 
the intended purpose; rather, the material presented here 
is to demonstrate how constraint damping and the gauge 
evolution equation (|16|l affects the evolution of the class 
of asymptotically flat, scalar field collapse binary black 
hole spacetimes considered here. 

The majority of results in this section will be from an 
equal mass head-on collision, though it will be shown that 
qualitatively similar conclusions apply to the orbital sce- 
narios described in Sec. The reason for looking at a 
head-on collision is that the simulation can be run in ax- 
isymmetry, and so consumes less of the limited computer 
resources that are available^. This makes it practical to 
do a more thorough survey of the effects of different con- 
straint damping and gauge evolution parameters. 

The parameters for the head-on collision, including ini- 
tial separation and simulation parameters (grid hierar- 
chies, dissipation, etc.) were chosen to be close to the 
parameters used in the 3D simulations. Units are scaled 
to Mq = 2mo, twice the initial mass of one of the black 
holes in the binary as measured by the area of its appar- 
ent horizon®. One common convention in the literature 
is to scale by the ADM mass of the spacetime; that is 
not used here because it includes the part of the scalar 
field that does not fall into the black holes, which could 
contribute as much as 10 — 15% to the total mass of the 
spacetime. The initial conditions for the scalar field in 
the the head-on collision simulation are as follows. The 
scalar field pulses are initially at rest (i.e. zero boost), 
and placed such that the initial coordinate separation of 
the two coordinate centers of the AH's that are first de- 
tected (at t « 2. 5 Mo) is 8.90Afo ^ the initial proper sep- 
aration measured along a coordinate line connecting the 
centers of the AH's from the surface of one AH to the next 
is IO.SA/q. With these initial conditions the black holes 
merge at i ~ 37Mo, and the energy emitted in gravita- 
tional waves is approximately O.OOIOA/q, calculated using 
(|31|l at a radius r = 50Mo. The "canonical" simulation 
relative to which others will be compared uses a value of 
the constraint damping parameter k = 1.80/Afo H14|) . and 
gauge evolution parameters = 0.50/Mq, £_2 = 5.4/Mo 
and T] = 1 H16|l . Factors of l/Afg have been inserted ac- 
cording to the dimensionality of the terms the constants 
multiply in the equations. There is no particular rea- 
son why these specific numbers were chosen, and, as will 



for example, a typical axisymmetric calculation takes on the or- 
der of an hour to several days, running on between 4 and 16 
Pentium IV class nodes of a Beowulf cluster. The 3D calcula- 
tions take on the order of several days to a month to complete, 
running on between 32 and 128 Pentium IV class nodes 
in the units were scaled by mo instead of 2mo 
the coordinate centers in unsealed units are at x = ±0.08, where 
X runs along the axis of symmetry. The tangent compactification 
function used maps x = ±1 to spatial infinity, thus x = ±0.08 is 
well within the linear range of this transformation |14| 



be demonstrated below, no fine-tuning of the parameters 
is necessary — the only requirement is that in magnitude 
the constants be of order unity relative to the smallest 
relevant length scale in the problem (Mq) to produce a 
noticeable effect. 



A. Convergence 

For convergence tests four characteristic grid resolu- 
tions are used, summerized in Table ^ below. In ax- 
isymmetry computational resources are available to go 
to higher resolution, though this has not been done to 
facilitate comparison with the 3D simulations, which to 
date have only been run with resolutions comparable to 
the three lower resolutions in Table |l| The first measure 
of the accuracy/convergence properties of the solution is 
shown in Fig. 2] which plots the sum of AH masses of 
black holes as a function of time. Assuming the AH is 
a good approximation to the event horizon, the sum of 
AH masses should approach a constant after scalar field 
accretion ends, and at late times after the merger. As 
seen in the figure, as resolution increases the conserva- 
tion of AH mass improves where it is expected to do so. 
Around the time of the merger there is a sharp spike-like 
feature in this function. This is a reflection of the robust- 
ness problems in the AH finder as discussed in Sec. lHl CI 
however, the amount of time that this anomalous behav- 
ior persists does seem to converge away. 

Fig. 12 shows TZh, an £2 norm of a residual of the Ein- 
stein equations in original formQ. Note that by moni- 
toring this particular residual rather than only the con- 
straints (|13|) , or the residual of the equations in GH form 
(|14|l . one has an additional check that errors have not 
been introduced in going from (O to ifTHl . and that the 
requirements for solutions of the GH form of the equa- 
tions to be solutions of the Einstein equations are sat- 
isfied. TZh is computed as follows. First, the residual 
T^abihj, k)h of all ten field equations 

Rab - 8^ {Tab - \9abT^ (34) 

is calculated using standard second order accurate finite 
difference approximations from a numerical solution ob- 
tained with a characteristic discretization scale /i, at a 
given grid location (i,j,k) (or in axisymmetry) . 

The residual at each point is normalized by the £2 norm of 
all second derivatives of all metric elements at the same 
point. This somewhat arbitrary normalization is sim- 
ply to give a convenient scale to plot the residual; the 
numerical value of the residual from one simulation is 
not particularly meaningful, rather it is the convergence 
to zero of the residual with increasing resolution that 
is a test of the correctness of the solution scheme. Af- 
ter computing the normalized TZab(,h 3,k)h, the infinity 
norm over the ten residuals is taken to define the resid- 
ual TZ{i,j, k)h of the grid location (i, j, k). The quantity 
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TZfi, as shown in Fig.|21 is then the £2 norm of TZ{i,j, k)h 
over the domain, though for computational convenience 
the points included in the norm were restricted to a uni- 
form distribution of size 129^ (129 x 65 in axisymmetry) 
that encompassed roughly bQM^ of the domain centered 
about the origin. The residual outside this region drops 
to zero quite rapidly. At early times (the first ~ 40Mo) 
the residual is dominated by the scalar field — i.e. it is 
largest in the vicinity of the outgoing waves of scalar 
field that did not fall into the black holes — which is the 
reason for the relatively large values of the residual then. 

For a finite difference numerical solution that is in the 
convergent regime one expects any quantity Qh (t) , calcu- 
lated from a solution obtained with a discretization scale 
h, to have a Richardson expansion of the form 



Qh{t) = Q{t) + eQi{t)h^ + eQ2{t)h^ 



(35) 



where Q{t) is the continuum value, eQi(i), eQ2(i)j •■• are 
a set of error functions that depend on Q{t) though not 
the resolution, and n is the order of convergence of the 
numerical technique^". With second order accurate dis- 
cretization n = 2. For the residual TZh{t), the continuum 
value TZit) — 0, and so ignoring higher order terms one 
has 



(36) 



Using residuals from simulations with two resolutions hi 
and /i2 one can eliminate the unknown error term from 
the above equation, and solve for n: 



n{hi,h2){t) 



\og{nhAt)/T^hAt)) 

log(/li//l2) 



(37) 



A measurement of n{hi, h2)(t) can be used as a con- 
vergence test: n{hi,h2) tending toward 2 as resolution 
increases implies that the assumed expansion in (|35|l is 
valid and the ignored higher order terms are small. If 
the continuum value Q{t) is not known, one can esti- 
mate the error in Qh2(t) calculated from a simulation 
with mesh spacing /12 if a second simulation with spac- 
ing hi is available, and assuming both simulations are in 
the convergent regime: 



Q{t) - QhAt) - eQ,it)h'^ 



h'l 



(38) 



Fig. |31 below shows n calculated for the sequence of 
residuals shown in Fig. [5] As the resolution of the pair 
of mesh spacings (/ii,/i2) increases, one expects the nu- 
merically calculated n to tend to 2; this trend is evident 
in the figure, especially at later times. At early times 



In an adaptive solution scheme one does not have a single dis- 
cretization scale h, nevertheless one would still expect a similar 
expansion with some characteristic scale h describing how well 
features of the solution are resolved 



"Resolution" 


wave-zone res. 


orbital-zone res. 


BH res. 


h 


1.7 Mo 


0.23Mo 


0.057Mo 


6/8 h 


l.SMo 


0.17Mo 


0.043Mo 


4/8 h 


0.85Mo 


0.12A/0 


O.O29M0 


3/8 h 


0.64Mo 


0.087Mo 


O.O2IM0 



TABLE I: The four sets of characteristic resolutions used in 
simulation results presented here, were each resolution is la- 
beled relative to the coarsest resolution h (only the three lower 
resolutions have been used for non head-on mergers). The 
grid is adaptive with a total of 8 levels of refinement, and the 
coordinate system is compactified. The wave zone is defined 
to be at r = 50Mo, the orbital zone within about r — IOMq 
and the black hole zone is within 2 — 3Mo of each apparent 
horizon. A CFL (Courant-Friedrichs-Lewy) factor of 0.15 was 
used in all cases. 
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FIG. 1: The normalized sum of total AH mass l|29^ as a func- 
tion of time, for the head-on collision simulations described in 
Sec. IIVI with constraint damping (compare Fig^J. This plot 
demonstrates convergence to a conserved mass before and af- 
ter the merger. The "spiky" behavior about the merger point 
(around t — 37 Mq) is due to AH finder problems as discussed 
in Sec HTTCl 



this behavior is not as apparent, though again due to 
the scalar field. Relatively speaking, the scalar field is 
under-resolved compared to the smallest scale of interest 
(the black holes): the smallest length scale in the scalar 
field is A H28|l . which was chosen to be roughly Mo/2, 
whereas the corresponding length scale for a black hole 
is its diameter, which initially is AniQ = 2Mo. For these 
simulations the scalar field is merely a vehicle to produce 
black holes, so that it is somewhat under-resolved is not 
a source of much concern. 
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FIG. 2: A norm of the residual of the Einstein equations 
@, calculated as discussed in Sec. IIV A1 for the head-on col- 
lision simulations described in Sec. II VI with constraint damp- 
ing (compare FigEl though note the different vertical and 
horizontal scales). This demonstrates the convergence of the 
solution — see also Fig. El 



FIG. 4: The normalized sum of total AH mass H29^ as a func- 
tion of time, for the head-on collision simulations described 
in Sec. II VI and without constraint damping; compare Fig. Q 
The curves end when the simulations crashed. 



B. Constraint Damping 




t/Mo 



This section contains a few results demonstrating the 
effectiveness of constraint damping. Fig. ^ is the equiv- 
alent of Fig. ^ showing the sum of apparent horizon 
masses as a function of time, though now the correspond- 
ing set of simulations have been run with the constraint 
damping parameter k = 0. All of the simulations without 
constraint damping "crashed" before the holes merged, 
which is why the curves end abruptly. Fig. |31 shows plots 
of the residual TZh for the simulations without constraint 
damping — compare to Fig. [3 though note the different 
scales used in the two plots. This clearly shows how effec- 
tive constraint damping is, though does not tell us how 
large k needs to be to start to work. To get an idea 
for what the required magnitude of k is, a set of 4/8/i 
resolution simulations were run with a range of different 
k's — Fig. shows the residual from these simulations. 
This demonstrates that k needs to be of order unity (rel- 
ative to 1/Mq), though no fine tuning is necessary, i.e. 
there are a range of values n > kq that are all essentially 
equivalent in damping rapid growth in the constraints. 



FIG. 3: The order of convergence 1371 calculated from the 
head-on collision simulations described in Sec. II VI (with con- 
straint damping), and using the residual data shown in Fig.|5| 
The discretization scheme is second order, and so one would 
expect n(/ii,/i2) to asymptote to 2 as the resolution is in- 
creased. This is evident in the figure, particularly after most 
of the scalar field has left the domain (around t « 40 — 50Mo). 



C. Source Function Evolution 

Fig. [7| demonstrates the effect of changing parameters 
in the source function evolution equation (|16|l for the 
head-on collision simulation; Table ^1 shows the values of 
the parameters for each set of curves plotted. To avoid 
clutter in the figure results are given from only 4 differ- 
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FIG. 5: A norm of the residual of the Einstein equations |(2J, 
calculated as discussed in Sec. If V Al for the head-on collision 
simulations described in Sec. \IV\ without constraint damping; 
compare Fig|21 though note the different vertical and horizon- 
tal scales. Again, convergence is evident, though the rapid 
growth of the residual with time prevents useful results from 
being obtained at modest resolution. 



1 — 1 — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 



1 - 




«=0 

K=0.18/Mo 

- K=0.31/Mo 

- K=0.45/Mo 
■ «=0.90/Mo 

K=l.B/Uo 

- /c = 3.6/Mo 







50 



100 

t/M, 



150 



200 



FIG. 6: Norms of the residual of the Einstein equations 
@, calculated as discussed in Sec. If V Al for a set of iden- 
tical "4/i/8" resolution head-on collision simulations, though 
with differing values of the constraint damping parameter k. 
The results shown in Fig. |5|were obtained with k — 1.8/Afo 
(though note the different vertical scale). 
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TABLE II: Values of the parameters in the gauge evolu- 
tion equation I16II used for the comparison simulation results 
shown in Fig. Q 



ent simulations, however these are quite representative 
of the range of behavior seen with this gauge evolution 
equation. The figure shows the lapse function a along the 
axis of symmetry at a few select times (at t = 0, a « 1 
and is identical in all cases). There are several features 
shown in the figure worth mentioning. First, the curve 
{D) demonstrates source function evolution without a 
damping term; this tends to cause oscillations about 1 in 
the lapse, and with sufficiently large values of the os- 
cillations grow to the point where the code crashes. This 
is the reason the (D) curve is not visible at late times. 

Second, in cases with a source term {A,C,D) a does not 
evolve to 1 in the strong-field regime even at late times, 
though is closer to 1 than harmonic gauge (B). For ex- 
ample, near merger at t — 29.9A/o the minimum value 
ttmin outside the excision surface for harmonic gauge is 
ttmin = 0.213 compared to 0.505 for case C, and at 
t = 200M anun = 0.458 (harmonic) versus 0.607 (C). 
These differences might not seem too significant, and at 
least for head-on collisions all cases except D work ade- 
quately. With non head-on collisions where there is sig- 
nificant orbital motion prior to merger, the small values 
that the lapse tends to with harmonic gauge seem to be 
problematic in that instabilities form before a common 
horizon is detected. There is not convincing evidence 
that this is coordinate problem as opposed to a numeri- 
cal issue, though given how "expensive" simulations are 
in 3D and that the non- harmonic gauge works, this issue 
has not yet been explored in any detail. 

Third, notice that for all gauge evolution schemes that 
are long-time stable {A, B, C), at late times the solution 
approaches a largely time independent form. This gauge 
condition therefore seems to be "symmetry-seeking" , at 
least for this class of spacetimes. 



D. Comparison with 3D evolutions 

Here it is demonstrated that some of the results just 
discussed for the 2D, axisymmetric head-on collisions ap- 
ply to the full 3D merger problems described in Sec.0 by 
comparing a select case. Fig. |S1 below shows the residual 
TZh from three 3D simulations with similar grid structure 
to the three lowest resolution 2D results shown in Fig.|21 
the sum of AH masses with time is shown in Fig. ^ and 
though there is no direct comparison to the 2D results 
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FIG. 7: The lapse function a along the axis of symmetry 
recorded at several times from a set of "4/i/8" head-on col- 
lision simulations as described in Sec. IIVI with different pa- 
rameters in the evolution equation 11611 used for Ht — see Ta- 
ble |n] Case A features the same values as the simulation 
results shown in Fig.'s ITISI case B is pure harmonic gauge, 
case C has a larger 77 coefficient than case A and so is more 
effective at keeping the lapse closer to 1, and case D is the 
same as A though with zero damping coefBcient ^2. Without 
the damping term the driving term, controlled by ^1, is too 
large in this example, and a overshoots 1 by a sufficiently 
large factor that the simulation crashes, which is why curve 
D is not present in the last two frames. 
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FIG. 8: Comparison plots of the norm of the residual of the 
Einstein equations (|5J , calculated as discussed in Sec. IIV Al 
between the head-on collision simulations described in Sec. IIVI 
(labeled 2D— see FigHJ and non head-on mergers (labeled 
3D) as discussed in Sec. [3 3/i/8 resolution simulations were 
not run in 3D due to lack of computational resources. Note 
also that the momentary increase in the residual near t = 
lOOMo — 150Afo in the 3D simulations is around the time of 
merger. 

Fig. ^1 shows the estimated Kerr spin parameter a after 
the merger. The parameters for the 3D simulations were 
chosen so that the equal mass merger occurs in roughly 
an orbit and a half for each resolution — see Fig. ^2 Note 
that the boost parameters differ slightly for each resolu- 
tion; this, as discussed in SecO is due to the sensitivity 
of the resultant orbit to the initial data parameters. Of 
course, in the continuum limit one would expect conver- 
gence to a definite orbit for a given boost parameter, or 
conversely, given a desired orbit the required boost pa- 
rameter should converge to a definite value. However, 
what Fig.|Slis meant to show is the behavior of the resid- 
ual TZh with time and as resolution varies, and so for 
a meaningful comparison simulations with similar grid 
structures and run-times were used. 



V. SCALAR FIELD COLLAPSE DRIVEN 
BINARY BLACK HOLE SPACETIMES 

In this section some results are presented from an ongo- 
ing study of scalar field collapse driven black hole merger 
simulations. The disadvantages of this approach to con- 
structing an orbit, from the point of view of simulating 
astrophysically relevant binary configurations, is that one 
does not have a priori control over the orbit that will re- 
sult, and perhaps more problematic is that in the regime 
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FIG. 9: The normalized sum of total AH mass 12911 as a func- 
tion of time for 3 merger simulations (see Fig llH , as described 
in Sec. |3 In each case the boost parameter was chosen so 
that roughly one and a half orbits is completed before merger, 
to better facilitate the comparison of the effect of resolution 
on solution accuracy (see the discussion in Sec.0. The late- 
time problems with the AH finder for the 6/i/8 case seems to 
set in when numerical error causes the angular momentum of 
the black hole to approach extremality-see Fig. IIUI 



where the evolution is started the radiation-reaction is 
sufficiently strong that it is not easy to deduce an astro- 
physical orbit that the binary might have evolved from. 
Therefore, these binary might be in a region of param- 
eter space that are not shared by astrophysical systems. 
However, from a theoretical perspective of trying to un- 
derstand the dynamics of binary black hole coalescence 
in general relativity, this scalar field collapse initial data 
is a straight-forward and simple way to probe a large and 
interesting region of parameter space. 

Here the focus will be on one particular family of ini- 
tial conditions: identical initial scalar field distributions 
{A = 0.3 and A = 0.827Mo) except that one is located 
at {x, y, z) — (4.45Mo, 0, 0) and given a boost with boost 
parameter v in the positive-y direction, while the other is 
located at {x,y,z) = (— 4.45Mo, 0, 0) and given a boost 
V in the negative-y direction — see (|28f) and the discus- 
sion afterwards. Approximately 85% of the energy in 
the scalar field falls into the black holes that form, while 
the remaining energy radiates away on roughly an orbital 
light-crossing time scale. 

This choice of initial data gives a single parameter v 
that can be varied. The black holes are initially suf- 
ficiently far apart — a coordinate (proper) separation of 
8.90Afo (IO.8M0) — that with a large enough boost pa- 
rameter one would expect the system to be unbound. 
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FIG. 10: Estimates of the angular momentum of the fi- 
nal black hole, corresponding to the simulations shown in 
Fig.'sF^and llll calculated using Iji-iO^ and the dynamical hori- 
zons framework in the top and bottom plots respectively. To 
facilitate the comparison the time has been shifted so that 
t = corresponds to the moment an encompassing AH is 
first detected. Note that with increased resolution the drift 
in a/A'I decreases at late times, as to be expected in a conver- 
gent solution. The 6/1/8 resolution case demonstrates what is 
seemingly a further source of robustness problems in the AH 
finder when the angular momentum approaches extremality 
(the spin- up is numerical-error driven here). 



At the other extreme, w = will give a head on collision. 
Thus there is a range of values < Vc for some v = Vc that 
will result in a merger, and it is this region of parameter 
space we want to explore. For the purposes of studying 
the gravitational radiation produced by mergers, and in 
particular comparing the orbital to plunge and ringdown 
part of the waves, it is useful to have parameters that re- 
sult in several orbits before the merger-see II 21 for a plot 
of the orbit from one such parameter. This regime is 
where these BBH simulations become quite interesting, 
in that there is apparently extreme sensitivity of the solu- 
tion as a function of the boost parameter. In particular, 
the dependence of the number of orbits n as a function of 
f , n{v) is quite reminiscent of certain aspects of critical 
phenomena in gravitational collapse in that there are ap- 
parently two distinct end states, merger of the black holes 
for V < V* versus separation for v > v*, and the closer 
V is to V* the more orbits before merger/separation that 
are observed — see Table iHll below for data from several 
simulations. Fig. ^|for a plot of the orbits from two of 
the 6/8/1 resolution runs that have been fine-tuned the 
most to date, and Fig. E| for graphs of the correspond- 
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ing gravitational waves emitted"'^^. This behavior is also 
quite reminiscent of "zoom-whirl" orbits seen in the tra- 
jectories of point-particle orbits about black holes j4l|. 
and might simply be the fully non-linear equivalent of 
the zoom- whirl phenomena. In fact, for zoom- whirl par- 
ticle orbits there is exponential sensitivity of the number 
of orbits iV„, completed in each whirl phase to the initial 
conditions : p—ps = e~°^™ , where p is the semi-latus rec- 
tum of the orbit, Ps is the separatrix of bound orbits, and 
a is some constant that depends (amongst other things) 
on the eccentricity of the orbit and spin of the central 
black hole^^. Of course, in the full solution where en- 
ergy is lost through gravitational wave emission can 
not be arbitrarily large. Furthermore, v* loosely defined 
above cannot be a "true" boundary between merger and 
separation, for one would still expect the orbits to be 
bound for v slightly greater than v* — the binary will just 
"whirl" out some distance before later merging. 

Note that the sensitivity to initial conditions compli- 
cates convergence testing (i.e. verification) of this inter- 
esting indication that zoom-whirl like behavior is present 
in equal mass binary merger spacetimes. This is because 
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FIG. 11: A plot of the orbits of the 3-D simulations used to 
compare with the head-on collision results (Fig. (HJ. What is 
shown for each case is the coordinate position of the center of 
the apparent horizon of one of the black holes. 



note however that no claims are being made that this is critical 
phenomena in any traditional sense of the word; rather, some of 
the terminology is useful in describing these tentative results, as 
is the bisection search strategy used in critical collapse studies 
to find the threshold parameter v* 



FIG. 12: A plot of the orbit of the v = 0.21906, 6^/8- 
resolution merger simulation. The spiral trajectories are the 
coordinate positions of the centers of the apparent horizons 
of the two black holes before merger, and the labeled curves 
show the coordinate shapes of the AH's, in the 2 = plane, 
at select times. 



sensitivity to initial conditions implies that the "bifur- 
cate" regime of parameter space must be found inde- 
pendently using several resolutions; i.e. one cannot do 
a "quick" low-resolution search to find v ^ v* , then 
perform a high resolution simulation to verify the low- 
resolution result, as each resolution will have a different 
V* . Furthermore, the numerical error that accumulates 
per orbit, as judged using AH mass conservation, is rel- 
atively large, in particular for the lower resolution sim- 
ulations, where the AH mass changes by a few percent 
per orbit. This is larger that the energy lost through 
gravitational wave emission, and though it does not ex- 
actly make sense to equate numerical error with physical 
energy, that these two "effects" are of the same order 
of magnitude suggests one cannot yet draw significant 
conclusions from these early results. 



VI. CONCLUSIONS 

In this paper further details were given of the gen- 
eralized harmonic evolution scheme introduced in [lj| 
and shown to be capable of simulating binary black 
hole coalescence Q. In particular, topics included were 
a demonstration of the effectiveness of constraint damp- 
ing in a fully non-linear setting, examples of the effect 
of source function evolution on the time slicing of the 
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h-resolution runs 



V 


n 


Pm/Mo 




ruf/Mo 


a/nif 


(E/Mo) 


0.21000 


1.3 






0.89 ±0.03 


0.75 ±0.05 


0.032 


0.21125 


1.4 






0.88 ±0.03 


0.74 ±0.05 


0.035 


0.21234 


2.3 






0.83 ±0.03 


0.73 ±0.05 


7 


0.21250 


2.7 


4.0 


3.6 






0.020 


0.21500 


1.5 


5.5 


4.6 






0.006 


0.22000 


1.0 


7.2 


5.8 






0.005 



6/8 h-resolution runs 



V 


n 


Pm/Mo 


dm /Mo 


■nif/Mo 


a/mj 


{E/Mo) 


0.20960 


0.9 









97 ±0.01 





65 ±0 


03 


0.028 


0.21750 


1.4 









92 ±0.01 





72 ±0 


03 


0.037 


0.21875 


2.0 









88 ±0.01 





70 ±0 


03 


0.046 


0.21906 


2.4 









86 ±0.01 





70 ±0 


03 


0.052 


0.219180 


2.8 









82 ±0.02 





70 ±0 


05 


0.063 


0.219200 


3.0 









80 ±0.02 





75 ±0 


05 


0.064 


0.219209 


3.3 









78 ± 0.02 





71 ±0 


05 


0.067 


0.219214 


3.7 









75 ± 0.02 





71 ±0 


05 


0.074 


0.2192175 


3.9 









74 ± 0.02 





70 ±0 


05 


0.076 


0.2192181 


4.3 









73 ± 0.02 





71 ±0 


05 


0.079 


0.2192188 


4.9 


3.2 


3.0 












0.058 


0.21938 


2.5 


4.8 


4.2 












0.019 


0.22000 


1.9 


5.3 


4.4 












0.014 



4/8 h-resolution runs 



V 


n 


Pm/Mo 


dm/Mo 


TO/ /Mo 


a/mj 


{E/Mo) 


0.21500 


1.4 






0.945 ±0.005 


0.71 ±0.02 


0.042 


0.22000 


2.1 


5.7 


4.8 






0.008 



TABLE III: Information extracted from the scalar field collapse driven binary simulations run to-date, described in Sec. [3 
Results from three different resolutions (see Table^J are shown (data from the h resolution v = 0.21234 run was lost, preventing 
the calculation of E in that case), v is the boost parameter for the initial scalar field pulses and n the number of orbits completed 
either before a merger, or at the moment the binaries reached the same coordinate separation as at t = 0. For binaries that did 
not merge, Pm is the minimum proper separation measured a.t t = const, along the coordinate line between, though exterior to 
the AH's, and dm is the minimum coordinate distance between the coordinate centers of the AH's. For binaries that merged, 
the final mass and angular momentum of the black holes, as estimated from AH properties, is listed. E is the total energy 
emitted in gravitational waves, calculated as described in Sec. lIIIj^ at a coordinate radius of 50Mo from the origin. Note that 
the quoted uncertainties in the final BH properties only include an estimated uncertainty from looking at plots such as Fig. |3] 
and 1101 To obtain an estimate of the numerical errors, one can use the Richardson expansion. For example by comparing the 
n = 1.4 results between the 6/8/i and 4/8/i runs, and using 13811 . the estimated errors in the v = 0.21500 4/8/i simulation are: 
nif = 0.95 ± 0.02Mo, a — 0.71 ± 0.01 and E — 0.042 ± .004. Once a larger set of simulations have been completed at higher 
resolution a more thorough analysis of immerical errors will be performed. 



spacetime, a discussion of certain technical aspects of 
the code and problems with the robustness of the ap- 
parent horizon finder, and some results from an ongo- 
ing study of scalar field collapse driven black hole bi- 
nary simulations. There are many outstanding issues 
that need to be explored, including the applicability of 



constraint damping to more general scenarios and alter- 
native evolution schemes, gaining more evidence for (or 
against) the zoom-whirl type behavior seen here, study- 
ing a broader range of binary black hole merger initial 
conditions, investigating a larger class of source func- 
tion evolution equations (in particular to include non- 
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FIG. 13: A plot of the orbits of two of the 6/i/8-resolution 
simulations that have been tuned closest to the bifurcate-like 
point in boost-parameter space (see Table 11111 . The merger 
case orbit ends when an encompassing AH is detected. 



trivial spatial source functions) , and beginning to extract 
some information from the non-perturbative regime of 
the merger to aid in gravitational wave detection efforts. 
The scalar field driven binaries are arguably not too use- 
ful in regards to this latter point, primarily because of 
the difficulty in mapping the resulting binary to (some 
approximation of) an astrophysical binary. Nevertheless, 
general features of the waveforms could be studied. Fur- 
thermore, what these simulations suggest is that to even- 
tually obtain accurate waveforms in the most interest- 
ing cases will require significantly improved accuracy in 
the simulations. Given that the present simulations al- 
ready utilize significant computer resources, rather than 
increase the resolution a more realistic near term solu- 
tion would be to incorporate higher-order finite differ- 
ence techniques or spectral methods (0, already use 
fourth order differencing in space and time, and use 
fourth order in space and second order in time). Modulo 
the accuracy issues, numerical relativity finally seems to 
have entered the era where it will begin uncover what 
will hopefully be the very rich and interesting landscape 
of black hole interactions. 

Acknowledgments: FP gratefully acknowledges re- 
search support from CIAR. The simulations described 
here were performed on the University of British 
Columbia's vnp4 cluster (supported by CFI and 
BCKDF), WestGrid machines (supported by CFI, 
ASRI and BCKDF), and Dell Lonestar cluster at the 
University of Texas in Austin. 



[1] F. Pretorius, "Evolution of Binary Black Hole Space- 
times", Phys. Rev. LeU.95, 121101 (2005) 

[2] M. Campanelh, CO. Lousto, P. Marronetti and Y. Zlo- 
chower, "Accurate Evolutions of Orbiting Black-Hole Bi- 
naries Without Excision", Phys. Rev. Lett.96, 111101 
(2006) 

[3] J. G. Baker, J. Centrella, D. Choi, M. Koppitz and J. 
van Meter, "Gravitational Wave Extraction from an In- 
spiraling Configuration of Merging Black Holes", Phys. 
Rev. Lett.96, 111102 (2006) 

[4] T. Nakamura, K. Oohara and Y. Kojima, Prog. Theor. 
Phys. Suppl. 90, 1 (1987) 

[5] M. Shibata and T. Nakamura, "Evolution of three- 
dimensional gravitational waves: Harmonic slicing case" , 
Phys. Rev. D52, 5428 (1995). 

[6] T. W. Baumgarte and S.L. Shapiro, "Numerical inte- 
gration of Einstein's field equations", Phys. Rev. D59, 
024007 (1999) 

[7] M. Campanelli, CO. Lousto and Y. Zlochower, "The 
Last Orbit of Binary Black Holes", Phys. Rev. D73, 
061501 (2006) 

[8] J. G. Baker, J. Centrella, D. Choi, M. Koppitz and J. 
van Meter, "Binary Black Hole Merger Dynamics and 
Waveforms", Phys. Rev. D73, 104002 (2006) 

[9] F. Herrmann, D. Shoemaker and P. Laguna, "Unequal- 
Mass Binary Black Hole Inspirals", gr-qc/0601026 (2006) 
[10] J. G. Baker, J. Centrella, D. Choi, M. Koppitz, 
J. van Meter and M. Coleman Miller, "Getting a kick 



out of numerical relativity" , 'astro-ph/0603204' (2006) 
[11] M. Campanelli, CO. Lousto and Y. Zlochower, "Gravi- 
tational radiation from spinning-black-hole binaries: The 
orbital hang up", gr-qc/0604012 (2006) 
[12] B. Brugmann, "Binary Black Hole Mergers in 3d Numer- 
ical Relativity," Int. J. Mod. Phys. D 8, 85 (1999) 
[13] B. Bruegmann, W. Tichy and N. Jansen "Numerical 
simulation of orbiting black holes", Phys. Rev. Lett. 92 
211101, (2004) 

[14] F. Pretorius, "Numerical Relativity Using a Generalized 
Harmonic Decomposition", Class. Quant. Grav. 22 425, 
(2005) 

[15] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen and 
O. Rinne, "A New Generalized Harmonic Evolution Sys- 
tem", gr-qc/0512093 (2005) 

[16] G. B. Cook and H. P. Pfeiffer, "Excision boundary condi- 
tions for black hole initial data," Phys. Rev. D 70, 104016 
(2004) 

[17] A. Buonanno, G. Cook and F. Pretorius, in preparation 

[18] D. Garfinkle, "Harmonic coordinate method for simulat- 
ing generic singularities", Phys. Rev. D65, 044029 (2002) 

[19] B. Szilagyi, B. G. Schmidt and J. Winicour, "Boundary 
conditions in linearized harmonic gravity," Phys. Rev. D 
65, 064015 (2002) 

[20] B. Szilagyi and J. Winicour, "Well-Posed Initial- 
Boundary Evolution in General Relativity", Phys. Rev. 
D68, 041501 (2003) 

[21] M. C. Babiuc, B. Szilagyi and J. Winicour, "Test- 



16 



ing numerical relativity with the shifted gauge wave," 
gr-qc/0511154 (2005) 

[22] M. C. Babiuc, B. Szilagyi and J. Winicour, "Har- 
monic Initial-Boundary Evolution in General Relativity," 
Phys.Rev. D73, 064017 (2006) 

[23] C. Gundlach, J. M. Martin-Garcia, G. Calabrese and 
I. Hinder, "Constraint damping in the Z4 formulation 
and harmonic gauge," Class. Quant. Grav. 22, 3767 
(2005). 

[24] O. Brodbeck, S. FrittelU, P. Hubner and O. A. Reula, 
"Einstein's equations with asymptotically stable con- 
straint propagation," J. Math. Phys. 40, 909 (1999) 

[25] C. Bona, T. Ledvinka, C. Palenzuela and M. Zacek, 
"General-covariant evolution formalism for Numerical 
Relativity", Phys. Rev. D67, 104005 (2003) 

[26] H. Friedrich, "Hyperbolic reductions for Einstein's equa- 
tions". Class. Quant. Grav. 13, 1451 (1996) 

[27] J. Balakrishna, G. Danes, E. Seidel, W. Suen, M. Tobias, 
E. Wang, "Coordinate Conditions and Their Implemen- 
tation in 3D Numerical Relativity" , Class. Quant. Grav. 
13, L135 (1996) 

[28] M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. 
PoUney, E. Seidel, R. Takahashi, "Gauge conditions for 
long-term numerical black hole evolutions without exci- 
sion" , Phys. Rev. D67, 084023 (2003) 

[29] L. Lindblom and M.A. Scheel, "Dynamical Gauge Condi- 
tions for the Einstein Evolution Equations", Phys. Rev. 
D67, 124005 (2003) 

[30] M. Alcubierre, A. Corichi, J. A. Gonzalez, D. Nunez, 



<0 



S - 



0.1 
0.05 


-0.05 

-0.1 
0.1 

0.05 



0.05 
-0.1 




-1 — I — I — r 



v=0.2192188 



I I I I 



I I I 



I I I I 



-I— I- 




_J I I I I L. 



[31 

[32; 

[33 

[34; 
[35; 

[36 
[37 
[38; 
[39 
[40 

[41 
[42; 
[43; 
[44; 

[45; 



B. Reimann and M. Salgado, "Generalized harmonic spa- 
tial coordinates and hyperbolic shift conditions," Phys. 
Rev. D 72, 124018 (2005) 

C. Bona, J. Carot and C. Palenzuela-Luque, "Almost- 
stationary motions and gauge conditions in general rela- 
tivity," Phys. Rev. D 72, 124010 (2005) 

C. Bona, L. Lehner and C. Palenzuela-Luque, "Geomet- 
rically motivated hyperbolic coordinate conditions for 
numerical relativity: Analysis, issues and implementa- 
tions," Phys. Rev. D 72, 104009 (2005) 

G. B. Cook, "Initial Data for Numerical Relativity", Liv- 
ing Rev.Rel. 3, 5 (2000) 

R. Arnowitt, S. Deser and C.W. Misner, in Gravitation: 
An Introduction to Current Research, ed. L. Witten, New 
York, Wiley (1962) 

H. Friedrich, "On the Hyperbolicity of Einstein's and 
Other Gauge Field Equations", Commun. Math. Phys 
100, 525 (1985) 

PAMR (Parallel Adaptive Mesh Refinement) and 
AMRD (Adaptive Mesh Refinement Driver) libraries 
(http: //laplace .physics .ubc . ca/Group/Sof tware .h tml j 
G. Calabrese, "Finite differencing second order systems 
describing black hole spacetimes", Phys. Rev. D 71, 
027501 (2005) 

A. Ashtekar and B. Krishnan "Isolated and dynamical 
horizons and their applications," Living Rev. Rel. 7, 10 
(2004) 

S.R. Brandt and E. Seidel, "The Evolution of Dis- 
torted Rotating Black Holes II: Dynamics and Analysis" , 
Phys.Rev. D52, 870 (1995) 

A. Nerozzi, C. Beetle, M. Bruni, L. M. Burko and D. PoU- 
ney, "Towards wave extraction in numerical relativity: 
The quasi-Kinnersley frame," Phys. Rev. D 72, 024014 
(2005) 

L. Smarr, in Sources of Gravitational Radiation, ed. L. 
Smarr, Seattle, Cambridge University Press (1978) 
J. Thornburg, "Finding Apparent Horizons in Numerical 
Relativity", Phys. Rev. D54, 4899 (1996) 
J. Thornburg, "Event and Apparent Horizon Finders for 
3 4-1 Numerical Relativity," gr-qc/0512169 (2005) 
C. Cutler, D. Kennefick and E. Poisson, "Gravita- 
tional radiation reaction for bound motion around a 
Schwarzschild black hole," Phys. Rev. D 50, 3816 (1994). 
K. Glampedakis and D. Kennefick, "Zoom and whirl: Ec- 
centric equatorial orbits around spinning black holes and 
their evolution under gravitational radiation reaction," 
Phys. Rev. D 66, 044002 (2002) 



100 200 

t/Mo 



300 



FIG. 14: Gravitational waveforms from the two simulations 
depicted in Fig. ^] What is shown here are the two po- 
larizations of one of the dominant spin weight —2 spherical 
harmonic components of ^'4 1321 1. calculated over a coordinate 
sphere a distance r = 50Afo from the origin, and normalized 
by rMo. 



